Mean-field Density Functional Theory of a Three-Phase Contact Line 
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A three-phase contact line in a three-phase fluid system is modeled by a mean-field density 
functional theory. We use a variational approach to find the Euler-Lagrange equations. Analytic 
solutions are obtained in the two-phase regions at large distances from the contact line. We employ 
a triangular grid and use a successive over-relaxation method to find numerical solutions in the 
entire domain for the special case of equal interfacial tensions for the two-phase interfaces. We 
use the Kerins-Boiteux formula to obtain a line tension associated with the contact line. This line 
tension turns out to be negative. We associate line adsorption with the change of line tension as 
the governing potentials change. 

PACS numbers: 05.70.Np, 65.40.gp, 68.05.-n, 68.35.Md 

Keywords: Line tension; line adsorption; three-phase contact line; diffuse interface model; mean-field density 
functional theory; phase-field model; successive over relaxation; triangular grid 



I. INTRODUCTION 

Studies of contact angle play an important role for the 
understanding of wetting phenomena in many systems, 
such as adhesives [l[ , liquid droplet spreading Q , and cell 
adhesion Q. Although contact angles can be measured, 
their theoretical computation can be complicated. They 
are affected by many factors, such as surface tension, line 
tension, temperature, composition of the system, and im- 
purities, especially surfactants. Here, we focus our at- 
tention on a three-phase fluid system (Fig. [TJ, where the 
three-phase contact line (briefly contact line) is the line 
where three interfaces and bulk phases meet. In this case, 
line tension is the excess grand potential per unit length 
of the contact line, which is a collective effect arising from 
inhomogeneities of intermolecular forces around the con- 
tact line, such as van der Waals, hydration, electrostatic, 
and steric forces (see [l|). The relevant forces can be 
short range or long range, the latter of which have 

been treated by the membrane method [9l4l5j or in terms 
of interacting surfaces (l6l - [20j . For a review see [2l|. In 
this paper, we deal only with short range forces so the 
problem can be formulated in terms of local densities. 

We model our system containing a three-phase contact 
line in the framework of general mean-field density func- 
tional theory by means of a diffuse interface model, where 
the imbalance of intermolecular forces is modeled by a 
potential function and a gradient energy of the chemical 
constituents. Thermodynamic-based functional theories 
incorporated with diffuse interfaces were first introduced 
by Lord Rayleigh [221 ], followed by many others [23l - |26j ] . 
They show good agreement with available experiments 
(see [H]]). For a comprehensive introduction of density 
functional methods to problems involving interfaces, see 
Rowlinson and Wido m [2 711 . Similar methods, known as 
phase- field models [23, [2{| , have been introduced to solve 
dynamical problems, such as moving boundary problems 
[30l - i32T |. Our model relates to a ternary solution (actu- 
ally a pseudo binary) and employs a different potential 
from that of Widom et al., Q, which is also a two-density 




FIG. 1. Geometry of a system with a three-phase contact line 
and three interfaces within a triangular prism. The contact 
line is a straight line perpendicular to the base and the cap 
of the prism, which are Neumann triangles, and each of the 
interfaces is also perpendicular to the lateral boundary of the 
prism. The regions divided by the three interfaces contain the 
three bulk phases, which are labeled by a, /9, and 7. 8 a , 8p, 
and 8-, are the equilibrium dihedral angles among the inter- 
faces, and (f> a , 4>p, <fi-y are the corresponding supplementary 
angles. L is the length of the contact line, and R a p, Rp-y, 
and Fiya are the distances from the contact line to the three 
lateral faces of the prism. We assume translational symmetry 
along L, so the problem is two-dimensional. 



model. Our potential is symmetric in the densities and is 
easy to relate analytically to measurable physical quan- 
tities in the far- field limit. 

We consider three bulk fluid phases in a multicompo- 
nent system. As illustrated in Fig. [1] the geometry of 
the system is a triangular prism. The three planer inter- 
faces a/3, /?7 and 7a meet at a three-phase contact line of 
length L and divide the system into three bulk phases a, 
/3, and 7, which subtend dihedral angles 8 a , Op, and 8 1 . 
Each of the interfaces is perpendicular to one of the lat- 
eral faces of the prism. The base and the cap of the prism 
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are Neumann triangles, which are perpendicular to the 
contact line and the three interfaces. The distances from 
the contact line to the lateral boundaries of the domain 
are R a /3, Rf)y, and R ja . LRij is the area of the interface 
ij . We treat this system in a regime where gravity is 
negligible. Due to the translational symmetry along L, 
the problem is effectively two-dimensional. Ultimately 
we consider the limit in which all Rij — > oo. 

Classically, the problem is usually treated by regarding 
the interfaces to be mathematical planes (zero thickness). 
Since the interfacial tension is the excess grand potential 
per unit area, the equilibrium angles can be obtained 
by requiring zero variation of the excess grand potential 
for an infinitesimal variation of the location of the three- 
phase contact line. The well-known result is 

sin^Q, svaOp sin# 7 

C/3 7 CTq 7 <7 a p 

where a a j3, cr,g 7 , and cr Q7 are the interfacial tensions. In 
this way, the interfacial tensions can be related to a Neu- 
mann triangle, whose three side are proportional the in- 
terfacial tensions and whose three angles are the supple- 
mentary angles of the three dihedral angles. For exam- 
ple, in Fig. [TJ <fip = 7r — 6 p. However, the classical model 
does not include the diffuse nature of the interface, nor 
possible complexity near the contact line. 

II. DENSITY FUNCTIONAL MODEL 

We are interested in a thermodynamically-based de- 
scription of a system which is inhomogeneous because of 
the interfaces and the three-phase contact line. We follow 
the thermodynamic methods of Gibbs, which amounts 
to choosing the grand canonical ensemble [13, P 228] in 
statistical mechanics. Thus, densities of chemical com- 
ponents as well as entropy density are allowed to vary, 
while the conjugate field variables are held fixed. Assum- 
ing that the grand potential of the entire system exists, 
the excess grand potential Vt xs due to the inhomogeneity 
of the system can be defined as 

n xs = n - n b , (2) 

where Q is the grand potential of the entire system and 
Of, is the sum of the grand potentials of the three bulk 
phases as if they shared the entire volume. Due to the 
homogeneity of the bulk phases, we have fif, = — pV, 
where p and V are the common pressure and the total 
volume of the bulk phases, respectively. Thus, 

n xs = n+ P v. (3) 

By convention [27|, ch 8], £l xa can be regarded as arising 
from two kinds of inhomogeneities, one associated with 
the contact line and the other associated with the inter- 
faces, i.e. 

^is = Lt + LR a po a p + LRp^ap^ + LR ia a iai (4) 



where r is the line tension, and <xy is the interfacial ten- 
sion of the interface ij far from the contact line. L and 
Rij are defined in Fig. [TJ According to this convention, 
r is defined as if the interfaces, with their far field values 
of aij, extend all the way to the triple line where they 
meet. The form (j4j) of excess grand potential is to be 
understood in the limit of all Rij —> oo. 

Following Rowlinson and Widom (2?J we assume that 
il xs can be expressed as the integral of a density of 
the excess grand potential, so 

V xs = lJ V>(x)dA (5) 

In mean field density functional theory, is assumed 
to be a functional of the number densities of the chemical 
components pi, i = 1, 2, • • ■ , c, for a c-component system, 
and Pc+i — s, the entropy density. Symbolically, 

^(x)=V[{Pi(x)} i=1 ,.., c+1 ], (6) 

which also depends on the set of conjugate field variables 
{p>i}i=l,... ,c+ij where pi,pL2, ■•• , Pc are chemical poten- 
tials and p c +i = T, the temperature, ipfe) is a function 
of densities and field variables plus a gradient energy cor- 
rection, 

V(x) = F ({ Pi (x)} i=li ... . c+1 ) + G ({V Pl (x)} 2=1 ... c+1 ) , 

(7) 

where F is a local density of the excess grand 
potential, an approximation sometimes called point- 
thermodynamics [27|, p 43] , and G is the density of gradi- 
ent energy, which is usually taken to be a linear function 
of the |Vpi(x)| 2 . The minimization of !1 IS is analogous to 
the minimization of the integral of a Lagrangian, whose 
role here is played by ^(x). Then, the terms in |Vpi(x)| 2 
play the role of kinetic energies and F plays the role of 
the negative of the potential energy. 

For a homogeneous bulk phase, there is no gradient 
energy and the excess grand potential density V( x ) 
/• ' ! ] where 

c 

F = ui + p = e - Ts P-iPi + V 

(8) 

= e — pi pi + p = (bulk phase). 

t=i 

Here u> is the uniform density of the grand potential; 
whereas, e, s, and pi are the densities of the internal 
energy, the entropy, chemical constituents that are uni- 
form in each bulk phase. Assuming the densities of the 
state variables for the inhomogeneous part of the system 
have a similar relation to those in the bulk phases, we 
approximate F for the entire system as 

c c+1 

F = e - \] jijpi - Ts + p = e-y j p,jpj + p, (9) 

i=l i=l 
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where e = e ({pi(x)}i=i,--- , c +i) is the non-convexified in- 
ternal energy as a function of the non- uniform c + 1 den- 
sities. Since p = p ({/ii}i=i,... , c +i) is the common pres- 
sure of the bulk phases, it only depends on the set {pi}- 
In general, e ({pi(x)}i = i r .. jC +i) is a non-convex function 
that has three potential wells and the bulk phases are 
given by a common tangent plane construction. Thus, 
F > because the terms — Ts — X)i=i PiPi + P represent 
the subtraction of the common tangent plane of the bulk 
phases from the non-convexified internal energy. There- 
fore, the three potential wells that correspond to the bulk 
phases are located at F = 0, where each is locally tan- 
gent to that plane. Note that e — Ts is the Hclmholtz 
free energy density, as for a bulk phase. By means of an 
approximation discussed by @, p 60] and M 

, we can 

reduce this model that depends on c+1 densities to an ap- 
proximate model that depends on only c densities, pi, p2, 
and p c . This amounts to assuming that de/ds = T, 
as it would for a bulk phase [35|. Thus, d(e — Ts)/ds = 0, 
so the form (|7|) of ip(x) is approximated by 

V(x) = F({ Pi (x)} i=1> ... , c ) + G ({VftCx)}^!,... , c ) , (10) 

where G, as a correction of F, is assumed to be only a 
function of the gradients of c densities as well. F also 
depends on the fields pi and T. 



A. Model for Uniform Molar Volume 

In this paper, we treat a ternary system under the con- 
straint of constant and uniform total molar volume. We 
obtain a tractable problem by introducing an explicit po- 
tential that is symmetric with respect to the three chem- 
ical components. 

For a ternary system, c = 3. Under the simplifying 
constraint of constant total molar volume, 



pi + p2 + p3 = p = constant. 



(11) 



With this constraint, the system we treat is actually 
a pseudo-binary system that can be described by two 
independent concentrations, say p\ and pi. Moreover, 
this constraint means that the conjugate thermodynamic 
variables of p\ and p2 are the chemical potential differ- 
ences Mi = pi—p3 and Mi = P2—P3, where the pi would 
correspond to a system with variable molar volume. In 
symmetric form, our potential is 

cv ^ p y^ (fit- Pa) 2 {Pi~ Pbf no , 
F{pi,p 2 , P3) = B > , (12 

i ' pi 

i=l ' 

where B is constant with the units of energy per unit vol- 
ume, and p a and pb are parameters (units of concentra- 
tion), which may depend on T and the Mi. By imposing 
the constraint p a + 2pb = p, we locate the three wells at 
symmetric positions. 

We introduce the notation Xi = p^J p as the mole 
fraction of chemical constituent i ranging from to 1, 




FIG. 2. Gibbs triangle. The 
summation of the distances 
from any point (Xi, X2, -X3) 
inside the triangle to the 
three sides of the triangle is 
equal to one, i.e. X\ +X2 + 
X 3 = 1. 




(a) a=2/3 



(b) a=l/6 



The three-fold symmetric potential with contours 

2/3, 



FIG. 3 

plotted on the base of the Gibbs triangle 
the three minima are located between vertices and the center. 

1/6, 



(a) For a 



(b) 



For a 



There is a local maximum at the center, 
the three minima are located between the mid edges and the 
center. 



a = p a j 'p, and b = pb/p = (1 — a)/2. The constraint (jTTJ) 
reduces to Xi + X2 + X3 = 1 . The potential (TT^j) can be 

expressed as 



F 
B 



= f(x l ,x 2 ,x 3 ) =J2( Xi ~ fl ) 2 (^ - 6 ) 2 



(13) 



i=l 



The function f{X\, X 2 , X 3 ) was originally introduced by 
Eldred (36[. In terms of independent variables, one has 
a two-variable function, 



f(X 1 ,X 2 ) = f(X 1 ,X 2 ,l-X 1 -X 2 



(14) 



The combination of (X1.X2, X3) can be illustrated by 
a point in the Gibbs triangle as shown in Fig. [21 The 
compositions of the bulk phases are (a, &, 6), (6, a, b), and 
(b,b,a). When a > 1/3, the three minima are located 
between vertices and the center of the Gibbs triangle, as 
illustrated in Fig. 3(a) The potential has a local max- 
imum at Xi = X 2 = X3 = 1/3. For a < 1/3, the 
three minima are rotated by 30° to positions between 
mid edges and the center of the Gibbs triangle, as illus- 
trated in Fig. |3(b)[ 

In the reduced approximation (1101) of ^(x), we assume 
that the gradient energy density is a linear function of 
the squares of the gradients of the mole fractions of each 
chemical component: 



|=.g(VX 1 ,VX 2 ,VX 3 ) = ^||VX l | 2 . 

i=i 



(15) 



where li are positive constants (with dimensions of 
length) associated with each chemical component. With 
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the constraint VX 3 = -VJCi - VX 2 , the form (JTSJ of 
gradient energy density turns into a two-variable func- 
tion, 

ff(VXi,VI 2 ) ee 0(VXi, VX 2 , -VXi - VX 2 ). (16) 

Thus, by inserting the explicit forms of potential (|13[) 
and gradient energy density ([T5j) into the form (fTOj) of 
V>( x )> the excess grand potential in our model §5§ be- 
comes 

Q xs =Bl[ [f{X 1 ,X 2 )+g{VX u VX 2 )]&A. (17) 

J A 



B. Euler-Lagrange Equations 

In equilibrium, we require 5fl xs = for infinitesimal 
variations of X\ and X 2 . To constrain total mole number 
of a finite system, we could add two Lagrange multipliers 
for Xx and X 2 to the integrand. However, we effectively 
work on an open system with infinite domain and fixed 
parameters fii and T, so particle conservation is not an 
issue and the Lagrange multipliers are effectively zero. 
From another point of view, the bulk phase is reached 
when the distance from the three-phase contact line to 
the boundary is large compared to the diffuse region of 
the contact line. This implies that the mole fractions 
should satisfy the boundary condition VXj-n = 0, where 
n is the unit outward normal to the physical domain. 
Thus, we obtain two coupled Euler-Lagrange equations 



<>J - {l\ + ^)V 2 Xi - £ 2 3 V 2 X 2 = 0, 
- {l\ + £j)V 2 X 2 - i\V 2 X x = 0. 



8X t 
df_ 
dX 2 



(18) 



C. Asymptotic Analysis in Far-Field 

In order to make a connection with the sharp inter- 
face limit of our mean-field density model, we consider a 
transition from phase f3 to phase a in the far-field regime, 
which is far from the three-phase contact line relative to 
the interfacial width. This is illustrated in Fig. 21 which 
corresponds to the transition from the minimum of one 
well to the minimum of another. Consistent with our 
potential function, X% — b is a constant in this region, 
which also satisfies the boundary condition WX 3 • n = 0. 
Therefore, the problem is essentially a one dimensional 
problem in a single variable, which we take to be X\. 

With X3 = b, we replace X 2 = 1 — b — X\, and sub- 
stitute V 2 Xi = d 2 Xi/ds 2 in the form of (fT?) of excess 
grand potential, where s is a coordinate perpendicular 
to the a/3-interface measured from /3 to a. The excess 
grand potential in the far field regime reduces to 



BLw 



H(X 1 



ai2 ( dXi 
2 V ds 



ds, (19) 



P (b,a,b) 




y (b,b,a) 



FIG. 4. Diagram of a transition between two bulk phases 
at a distance far from the three-phase contact line. In a far 
field limit, the transition from bulk phase j3 to bulk phase 
a is represented by the transition between two wells from 
minimum (6, a, b) to minimum (a, b, b) in the Gibbs triangle. 
This corresponds to the transition along a one-dimensional 
coordinate s perpendicular to the interface in spatial space. 
Here, w is the width of an area at a distance far from the 
contact line. 



where w is the width of an area in the far field regime as 
indicated in Fig. H H(X t ) = 2(X 1 - a) 2 (X 1 - b) 2 , and 
«i2 ee i\ + t\. The limits of integration are effectively 
from —00 to 00. 

In equilibrium, we require SQ XS = for an infinitesimal 
variation of X\ and obtain the Euler-Lagrange equation 
in the far- field limit: 



H'(X 1 )=a 12 



d 2 X x 
ds 2 



(20) 



Then, we multiply by dX\/ds and integrate to obtain 
H(X 1 ) = ^(^)\ (21) 



ds 



where the integration constant is zero because H(a) — 
H{b) — and the slope dX\/ds is zero for X\ = a and 
X x = b. 

By solving (|2 1 [) , we obtain the far-field solution for X\ 
at the a/3-interface, 



tanh 



(22) 



where we choose s = as X\ = (a + b)/2 and define the 
interfacial width parameter of the a/3-interface as 



(23) 



Of course, X 2 = 1 — b — X\, These analytic solutions 
were originally found by Eldred [HI . Analytical far-field 
solutions for density profiles and interfacial tensions for 
a symmetric three-phase contact line but a different po- 
tential were obtained by Szleifer and Widom Q. Note 
that when a = b = 1/3, the interfacial widths diverge. 
The contact line and the three interfaces vanish. In this 
case, the three chemical constituents have mole fractions 
of 1 /3 distributed uniformly over the entire system. 
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s (distance from a/3 interface) 
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FIG. 5. Asymptotic far-field solutions for the mole fractions 
Xi at the a/3 interfaces, for a = 2/3. s is the distance from /3 
to a perpendicular to the a/3 interface 



As illustrated in Fig. [5j when s is negative infinity, we 
have X\ — b and X2 = a, which indicates the /3 bulk 
phase. In contrast, when s is positive infinity, we obtain 
X\ = a and X2 = b, which refers to the a bulk phase. 
Similarly, we can apply the same analysis for the other 
two interfaces and obtain solutions for X\, X%, and X3 
in the far- held limit. 

The definition of interfacial tension is the excess grand 
potential per unit area of interface. Thus, by inserting 
the relation (|2"T1) that connects potential and gradient 
density into the excess grand potential (p~9|) . the interfa- 
cial tension of the a/3-interface can be expressed as 



Cq/3 




(24) 



BjaH I y/2H(X 1 )dX 1 . 

After integration, the interfacial tension of the a(3- 
intcrfacc is found to be 



|a-6| 3 n 

fa/3 = 7, By/ Oti2 



\3a-l\ 
24~ 



BJef + el (25) 



The interfacial tensions of the /37-interface and the 7a- 
interface can be obtained similarly. Consistent with the 
classical relation of equilibrium angles ([1}, the equilib- 
rium angles in our model obey 



sin ( 



VW+tl VW+t 



(26) 



III. NUMERICAL ANALYSIS FOR 
SYMMETRIC THREE-PHASE CONTACT LINE 

Due to the nonlinearity of the Euler-Lagrange equa- 
tions (1181) , one cannot obtain an analytic solution for the 
entire domain containing the three-phase contact line. 



We consider a simplified symmetric contact line centered 
in an equilateral triangular prism. 

Let £1 — £2 — £3 = £, where £ is a characteristic length. 
For convenience, we define the dimensionless coordinate 



r/£, and V 



J V . Thus, the dimensionless form 



of the excess grand potential (JT7J) is given by 



it 



BLt 2 



[f(u,v)+g(V / u,Vv)}dA\ (27) 



where A' = A/£ 2 , and, for convenience of writing, we 
define u = Xi and v — X2', then f(u,v) is the potential 
(H3J), and g{W'u, V'u) is the symmetric version of the 
gradient energy density ([T6]) . but in the form 

g(Vu, Vv) ee |V'u| 2 + \Vv\ 2 + V'u • V'v. (28) 

Similarly, the Euler-Lagrange equations (|18p can be di- 
agonalized and take the dimensionless form 



V /2 ATi 



2 5/ 18/ 



3dXi 3dX 2 



= 0, 



(29) 



3 dXi 3 dX 2 
The dimensionless interfacial width parameter is 

6> nt = V2/\a-b\. (30) 

Because of the three-fold symmetry of our system for a 
symmetric three-phase contact line, and the fact that the 
Laplacian operator is well-behaved on a triangular grid, 
we employ an equilateral triangular grid to resolve our 
special geometry. The computational domain is chosen 
as an equilateral triangle with physical dimension large 
compared to the dimensionless interfacial width pop , and 
the grid points are determined by filling out smaller tri- 
angles with non-dimensional length d! as shown in Fig. [5] 
There are N grid points on each domain edge, which is 
of length H' — (N — \)d! and perpendicular to one of the 
interfaces. The distance from the contact line to each 
of the outer edges is R' = (N - l)d'/(2 v / 3). For conve- 
nience, we make a special choice of grid points to allow 
the grid points to lie at important points of our system, 
such as the center of the contact line and the transition 
points of the far-field interfaces. To do this, the number 
of grid points on each edge is chosen to be N — 6m + 1 , 
where m is an integer. Therefore, the dimensionless size 
of each edge of the outer triangle is H' = 6md' and the 
dimensionless distance from the contact line to each edge 
of the outer triangle is R' = y/Smd'. 

The asymptotic far-field solution approaches an effec- 
tively one-dimensional two-phase problem. The interfa- 
cial width is small compared to the distance from the 
three-phase contact line. This setup makes the far-field 
solutions easy to apply at the boundary of the domain. 
Also, the corner regions of the large triangular domain 
approache the bulk phases, where our potential vanishes. 
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ya-interface 




/?y-interface 



FIG. 6. A small triangular grid. The triangular domain is 
filled out by small triangles with dimensionless edge length 
d! . For convenience, each edge of the outer triangular domain 
is perpendicular to an interface. Also, by specific choice, the 
geometric center and the interfacial centers of every edge are 
on grid points. The triangular grid follows the rule that the 
number of grid points on each edge is N — 6m + 1, where m 
is an integer. Here we take m = 1 for illustrative purposes 
only. 



A. Consistent Discretization 

To discretize the Euler-Lagrange equations (|2Q[) . for a 
symmetric three-phase contact line based on the triangu- 
lar grid in Fig. [51 we employ a variation of the discrete 
form of the excess grand potential to avoid inconsistent 
discretization of the potential of / and the gradient en- 
ergy g. We approximate u and v as planer functions in 
the region of each small triangle of the triangular grid. 
Then, the value of u and v at the central point n is de- 
fined as u n = J2 m ev(") u ™/3 and Y^meVM v ™/3, where 
is the set of vertices of the small triangle denoted 
by its center point n (see Fig. 7(a) ). The discrete form of 
the dimensionless excess grand potential (fT7|) is approx- 
imated by evaluating the integrand at the central point 
of each small triangle, 



nECP 



(/„ + g n )A 



(31) 



where CP is the set of the central points of the small 
triangles over the entire triangular grid in Fig. [5J A is 
the area of each small triangle; f n = f{u n ,v n ) and g n = 
g{(y'u)m(y'v) n ). Since u and v are approximated by 
planer functions, we obtain 



fjn 



2 



E 

< 3 -v k f + (uj 



[(« 



Uk) 



(32) 



Uk)(Vj -v k )] 



where PF" is the set of pairs of the vertices 



(Fig. 7(a) I of the small triangle n with edge d! . 

At equilibrium, we require 5fl' xs = for the discrete 
form (f3"Tj) of the excess grand potential. From the chain 




FIG. 7. | (a) | A small triangle in the grid of the physical 
domain. Each equilateral small triangle has edge length d! 
and is denoted by its center point n. V^ n ' = {/ii, /12, ^3} is 



the set of vertices of the small triangle denoted by n. (b) 
The nearest neighbors and the nearest center points for a 
site i in a triangular grid. NCP^ 1 ' = {ni, ri2, n-j, 124, ns, ng} 
is the set of the nearest center points for the site i, and 
NN^ = {mi, mi, m3, 1714, m 5 ,me} is the set of nearest neigh- 
bors for the site i. 



rule, this is equivalent to the vanishing of the sum of 
the variations of all of the unknowns (ui,Vi) for each 
internal site i of the triangular grid. According to the 
approximations (|32j) of the gradient energy density and 
requiring the coefficients of Sui and Svi to vanish, we 
obtain the discrete Euler-Lagrange equations for each site 
i (each internal grid point of the triangular grid) , 



(V 2 ^ - - 




1 



NCPV 



NCPV 



< 9v J NCPV 

2 far 



= 0, 



3 V dv 



(33) 



= 0, 



NCPW 



where ( df/du) NCpW = E„ e wcpw (df/du) n /6 and 

(df/dv) NCp(l} = E„eArcp(')( 9 // 9u )n/ 6 are the ad- 
ages of df/du and df/dv over t he six nearest center 
points of each site i, NCP {€) (Fig. |7(b)[ ). The approxi- 
mate Laplacian operators according to second order Tay- 
lor's series expansions are (V' 2 u)i = 4(u NN w — Ui)/d' 2 
and {V 2 v)i ee 4(tJ JVAr(s) - Vi)/d /2 , where u NN (i) = 

Emewjvw ""i/ 6 and f jvw<') = Emewjv(-) u ™/ 6 are the 
averages of u and v over the six nearest neighbors of each 

site i, NN® (Fig. [7(b)]) . 

Note that the discrete Euler-Lagrange equations (|3"3")l 
are similar to the analytic form ([2U| , except df / dX\ and 
dj ' jdX-i are replaced by the average values over the six 
nearest central points. After we apply the asymptotic far- 
field solutions as the boundary conditions of the system 
of algebraic equations for the triangular grid, there are 
( N — 2) (N — 3) algebraic equations for the whole domain. 



B. Successive Over-relaxation Method 

To solve the system of coupled algebraic equations, 
we apply the method of successive over-relaxation (SOR) 
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(371 138| . We define error equations for the diagonalized 
form of the discrete Euler-Lagrange equations (|3"3")l : 



(r u )i = (V /2 u), 
{r v )i = (V' 2 u)i - ~ 



A 


2 




3 






A 




(m 


3 




\duj 



( 2l 

dv 



NCPW 



NCPW \ OV / JVCPW. 

(34) 



where < A < 1 is an adjustable parameter used to 
implement our numerical technique, and {r u )i and (r„)j 
are the residues that we try to make as small as practical. 

In the form (jMJ) of error equations, (df /du) NCp{i) and 
(df/dv ) NCP (i) are polynomials of the mole fractions. To 
avoid the complexity of numerical calculation due to the 
nonlinearity of these terms, at the beginning, A is set 
to be zero. After solving this simplest version of the 
equation by SOR, we apply that solution as the initial 
values of SOR for new equations in which A is increased 
by a small fraction of 1, and solve the equations again. 
Then, we gradually enlarge A and repeat this procedure 
until A reaches one. 

In the updating process of SOR, we first input guessed 
numbers of ui and Vi as initial values into error equations 
(1341) for every site in the grid. Then, we update Ui and 
Vi for each site by 



old 



d' 2 

-q—{rv)i, 



(35) 



where q = 1.86 39]. We repeat this procedure until the 
values of Ui and Vi at every site converge. 

To check convergence, we study the norm of errors after 
every iteration. The norms are defined as 



M = »^and 



N tnt 



(36) 



Then the convergence criteria can be defined as \\r u \\ and 
||r„|| are simultaneously smaller than e, where e is a small 
number. Alternatively, this means ||u™ et0 — and 



.old I 



are simultaneously smaller than qd' 2 e/4. 
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(a) Contour plot of the numerical solution of X\ 
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(b) Profiles of X\ along x' = 



FIG. 8. (a) Contour plot of the numerical solution of X\ for 
a = 2/3, m = 180, and d! = 0.1. Here, N = 1081 and b = 
1/6. x = x/l and y = y/£ are the dimensionless Cartesian 
coordinates. The domain edge is H' = 108, which is large 
compared to the interfacial width parameter S' int = 2v / 2. The 
contours are evenly spaced from 10% to 90% of a — b. Note 
that the interfacial width defined by the difference between 



10% and 90% at the boundary, S[ 



10%-9(>%.b> 



is around 2. 20(5' „ 



whereas the interfacial width at the three-phase contact line, 
5 10 %_9o%,t) is around 2.6l5' int . |(b)| Pr ofiles of Xi along x' — 
(along the central vertical line of Fig. |8(a) I and the boundary 
shared with a and j3 phases of the numerical solution of Xi 
for a — 2/3, m = 180, and d' = 0.1. The diffuse region of the 
profile along x' — is slightly widened and shifted compared 
to the profile along the boundary shared with a and /3 phases. 



Contours and Profiles 



Here, we present a numerical solution obtained from 
SOR. The numerical input is a = 2/3, d' = 0.1, and 
m = 180. The error tolerance, e, is 10 -8 . Then the 
domain edge H 1 = 108 is large compared to the inter- 
facial width parameter S' int = 2\/2. Solutions for X 2 
and X$ are just the rotation of the solution of X\ by 
120° and 240 



The solution for Xi in Fig 
veals the nature of diffuse interfaces. There is bending 
and slight widening of the diffuse region for X\ near the 



three-phase contact line, which is quantified by the in- 
terfacial width defined from 10% to 90% isoconcentra- 
tion lines. The width at the boundary, 5'-, 
about 6.21 fa 2.20(51 
contact line, S' w% _ go9i , 



10%-90%,fc' 1S 

int . For comparison, the width at the 
- 7.37 pa 2.61^' nt , is about 20% 
is small compared to the distance from 
re- the outer domain boundary to the contact line along any 



larger. 5', 



10%-90%,t 



interface, which is R' ps 31.2 ps \l.Q5' int . Also, the pro- 
file at the contact line shifts its center compared to the 
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one at the boundary, as shown in Fig. |8(b)| Close to 
the boundary, the nearly parallel isoconcentration lines 
along the interfaces show that our domain size is close 
to the asymptotic regime, consistent with our intended 
boundary condition. 



IV. LINE TENSION 

A. Density Functional Model for Line Tension 

The numerical results in Sec. IIIII reveal the fact that 
the actual interfacial width increases slightly while ap- 
proaching the three-phase contact line. This result is 
different than that which would be obtained by extrapo- 
lation of the far-field solution. In this section, we study 
the line tension which is the excess energy per unit length 
associated with the three-phase contact line. By conven- 
tion, the line tension is defined in the form (j4|) of the 
excess grand potential. For a symmetric contact line, we 



let R = R a f 



R 



and a = a a p — 17/3^ 



'^fa- 



Ill terms of the dimensionless grand potential (|27[) with 
R' = R/l, the dimensionless line tension is given by 



T ' 55 BP = ~ 3i?V 



(37) 



where a' is the dimensionless form of the interfacial ten- 
sion in the far-field limit (f2l)j) for a symmetric contact 
line, i.e. 



a - b\ s 



-(- 

V2\2 



(38) 



From the symmetry of ip{x) for a symmetric contact line, 
the solutions for the three mole fractions have the prop- 
erty Xi(r,ff) = X 2 (r,9 - ^) = X 3 (r,8 + 2jL), which 
means each of them are given by only a rotation of 
2ir/3 or — 27r/3 from the others. Thus, the integration 
in the form (|T7|) of excess grand potential can be di- 
vided into three equal parts. By applying the boundary 
condition V X ■ h = 0, we find that / \ |V'Ai| 2 dA' = 
— J ^XiS7' 2 XidA' . Thus, the dimensionless line tension 
is given by 



t' =3 



(Xx - af{X x - bf - \x l V 2 X 1 



V2\a-b\ 3 R'. 



dA' 



(39) 



We find, however, that evaluation of the form (|39|) of 
dimensionless line tension is sensitive to the choice of 
boundary condition, which may result from the incon- 
sistency between the numerical evaluation of the excess 
grand potential and the analytic interfacial tension. In- 
stead, we use a formula for line tension derived by Kerins 
and Boiteux , which transfers the second term of the 
form (|39p of line tension into a surface integral and com- 
bines it with the first term. In this integral form, the in- 
tegrands will vanish at distances far from the three-phase 



contact line, which means it is insensitive to domain size 
for a sufficiently large domain. According to the Kerins- 
Boiteux formula, the dimensionless line tension in given 
by 

t'= [ [-f(u,v)+g(Vu,V'v)]dA'. (40) 

J A 



Numerically, we can discretize the integral in Eq. (|40 
by employing a triangular grid as in Fig. [51 so 



A, 



(41) 



n£CP 



where /„ and g n are defined in the discrete form pip of 
the excess grand potential and the approximation (|32[) 
of the gradient energy density. Then we can utilize the 
numerical method developed in Sec. IIIII to obtain the nu- 
merical value of dimensionless line tension in our model. 



B. Evaluation of Line Tension 

We perform a numerical evaluation of the integrand 
of the discrete form (|4Tj) of the Kerins-Boiteux formula. 
Fig.|H]shows a contour plot of the integrand on a logarith- 
mic scale. A similar plot on a normal scale can be found 
in Taylor and Widom (35j | . The integrand decays expo- 
nentially for the most of the domain. The major contri- 
bution of the integrand is approximately within the range 
from 10~ 2 to 10 -5 at a core region centered at the three- 
phase contact line with dimension of 2 to 3 times S' int . 
The minor contribution, which is considered to be from 
10~ 5 to 10~ 8 , is distributed outside the core region and 
along the three interfaces with a width of about 1.5S' int . 
The integrand in the rest of the domain is less than 10 
and is essentially negligible compared to the one close 
to contact line and along the three interfaces. Theoreti- 
cally, the potential and gradient energy density are zero 
within bulk phases and — / + g — > in the interfaces far 
from the core, so the small but non-vanishing values of 
— f + g along the interfaces results from numerical er- 
rors. Also, we can see that the contours of the integrand 
begin to bend at the far-field boundary of the domain, 
which may relate to the errors associated with apply- 
ing the far-field solution as the boundary condition for 
a finite domain. Note that the contours along the three 
interfaces are nearly parallel except close to the contact 
line and the boundary. This suggests that the numerical 
evaluation of the Kerins-Boiteux integral over these ar- 
eas leads to an error that is approximately proportional 
to R' . Our numerical results also show that when d! is 
smaller, the distribution of the integrand within the core 
region is sharper, with a slightly larger maximum value, 
and decays faster, which means that the integrands along 
the three interfaces and the boundary decrease when d! 
becomes smaller. So, we assume that the numerical er- 
ror of the evaluation of the Kerins-Boiteux formula is 
proportional to R' and depends on dl . 
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FIG. 9. A contour plot on a logarithmic scale at the integrand 
of the Kerins-Boiteux formula over a triangular domain with 
a threshold of 10 -8 . The parameters for this calculation are 
m = 180, d! = 0.1, and a = 2/3. The major contribution of 
the integrand with values from 10~ 2 to 10~ 5 is confined in a 
core region with a dimension of 2 to 3 interfacial widths, S' int , 
near the three-phase contact line. The minor contribution, 
which has values ranging from 10~ 5 to 10 -8 , is outside the 
core region and along the three interfaces with width of about 
l.58' int . This shows that the integrands within the bulk phases 
are significantly smaller compared to the core region and along 
the three interfaces. Also, the nonzero contours along the 
three interfaces are nearly parallel near the outer boundary. 



To test this, we use the Kerins-Boiteux formula to cal- 
culate values of line tension, r'. As shown in Fig. [TU1 
the t' value is nearly proportional to R' for each d! . We 
take grid spacings, dl — 0.05, 0.1, 0.2, and 0.4, and do- 
main sizes, R 1 w 31.2, 41.2, 52.0, and 62.4, which are 
relatively large compared to the size of the three-phase 
contact line, roughly 7.37 « 2.&\8' int . By linear extrapo- 
lation from the results in Fig. [TOl we find that the values 
of t' for different d' roughly meet at R' = 0. Moreover, 
from Fig. 111! we find that the dominant numerical error 
of t' comes from a d' 2 term for fixed R' values because 
the calculated t' is almost linear in d' 2 . 

On the basis of Fig. QJJ and Fig. [Til we assume that 
the numerical value of r' is a function of R' and d' of the 
form 

t'(R', d') ~ f' + Cl d' 2 + c 2 R'd' 2 + h(d' 2 ,R'), (42) 

where c\ and ci are constants, f ' is the line tension nearly 
invariant of the grid spacing d! and distance R! , and 
h(d' 2 ,R') represents terms of higher order than d' 2 and 
R' . From the expression (|42|) of the numerical line tension 
r'(i?', d'), we can eliminate approximately the numerical 
error which depends on R' by linear extrapolation of r' 
from various R' toward R 1 = as in Fig. QJJ and obtain 
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FIG. 10. Line tension r as a function of R for various d' 
values calculated by Kerins-Boiteux formula with a = 2/3. 
r' is calculated for d' = 0.05, 0.1, 0.2, and 0.4 and R « 
31.2, 41.2, 52.0, and 62.4. The dashed lines are the linear 
extrapolations of the line tensions from various R toward 
R = for each d' . 
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FIG. 11. Line tension as a function of d' 2 at various R values 
calculated by the Kerins-Boiteux formula with a = 2/3. r' is 
calculated for d' = 0.05, 0.1, 0.2, and 0.4 and R w 31.2, 41.2, 
52.0, and 62.4. 



a correction of t' at R! = , which is given by 

r' {d') ee t'{R' = 0, d 1 ) ~ r' + ad' 2 + h (d' 2 ), (43) 

where ho(d' 2 ) represents terms of higher order than d' 2 . 
The extrapolated results of Tq (|4"3")) are listed in Table UJ 
In addition, we can refine our result at R' = by 
using Richardson's extrapolation [4lll42T |. in which results 
for two successive d' values are used to eliminate the d' 2 
term. The first level of correction is defined as 



r[{d') 



Ar' Q {d')-r' Q {2d') 



f' + 



/ii(d' ; 



(44) 



where hi(d' 2 ) represents the terms of higher order than 
d' 2 . From the calculation of r[ in Table HI the line tension 
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TABLE I. Refinement of line tension r' based on the Kerins- 
Boiteux formula. Tq is the extrapolated line tension at R' — 
for various d! , where the dominant term is d' 2 . t[ is the first 
level correction of Tq by eliminating the d' 2 term. 





d'=0.05 


rf'=0.1 


cf =0.2 


d'=0A 




-0.07216833 


-0.07216553 


-0.07215562 


-0.07211596 


/ 

T\ 


-0.07216926 


-0.07216883 


-0.07216884 


N/A 



t' for a = 2/3 is approximated by 

f' ~ -0.072169, (45) 

where the uncertainty is in the hnal digit. It is well known 
both theoretically [13, and experimentally [l|, [2l[ that 
line tensions, unlike interfacial tensions, can be either 
positive or negative. Physically, a negative line tension 
means, for example, that the line of intersection of a ses- 
sile drop with a substrate would tend to expand P, 13] , 
but is ultimately limited by positive interfacial tensions. 

C. Scaling of the Density Functional Model for 
Line Tension 



In terms of independent variables, we define the two- 
variable functions 

j(Yi ,Y 2 ) = f(Y 1 ,Y 2 ,1-Y 1 -Y 2 ) (50) 

and 

3(VFi,VF2) = 5(VFi,VF2,-VFi- VY 2 ) (51) 

Similarly, we obtain a new expression of the dimension- 
less line tension (|4"0"|) in the form of the Kerins-Boiteux 
formula |40| |: 

f = TZE = I {-KYx,Y 2 ) + ~g{VY 1 ,VY 2 ))<iA. (52) 

Bt £ J A 

Note that the integral of fl xs and f are both independent 
of a and dimensionless. 

Based on the numerical methods presented in Sec. IIII1 
we can compute f and refine the result by Richardson's 
extrapolation as in Table [TTJ The refined f is 

f ~ -0.28868, (53) 

where the uncertainty is in the last digit. 



In our original way of scaling, we factored out Bl 2 
from the excess grand potential and also from the Kerins- 
Boiteux formula for line tension. Then, we calculated the 
integral in a dimensionless domain. However, our poten- 
tial is parametrized by a, which means that solutions 
of the Euler-Lagrange equations and the calculation of 
t' depend on a. Here, to elucidate the a-dependence of 
our model, we study the problem in a new framework by 
defining the following new scaled variables, 



= Xj-b = 2X, -{I -a) 
a — b 3a — 1 



(46) 



where X^=i Y = 1 and a ^ 1/3. In our model, the value 
of Xi is limited from 6 to a, so lj varies from to 1. In 
terms of the new scaling variables, the new dimensionless 
form of the excess grand potential (|17|) for a symmetric 
three-phase contact line is given by 



it i 



fi, 



BLl 2 



if(Yi,Y 2 ) + s(VYi, VY 2 ))dA, (47) 



where we define the following new scaled constants B = 
B{a - 6) 4 , P = £ 2 /(a - b) 2 , V = IV, A = A/£ 2 . The 
scaled potential and gradient energy density are 



and 



f(Y 1 ,Y 2 ,Y 3 ) = J2(Y) 2 (Y-l) 



Kvy 1 ,vy 2 ,vy 3 ) = ^|vy i | 



(48) 



(49) 



TABLE II. Refinement of the scaling line tension f based on 
Kerins-Boiteux formula, fo is the extrapolated line tension at 
R — for various d values, where the dominant term is d 2 . fi 
is the first level correction of fo by eliminating the d 2 term. 











d= 0.05 


d = 0.1 


d = 0.2 


d = 0.4 


fo -0.28866211 


-0.28862233 


-0.28846449 


-0.28781290 


fi -0.28867521 


-0.28867560 


-0.28868169 


N/A 



According to the new scaled constants, we find that 
the dimensionless line tension r' (140|) and its new scaled 
expression f (|52p obey the following relation 



T 

BP 



(54) 



which shows that r' is proportional to (a — 1/3) 2 as indi- 
cated in Fig. [T2l and r' is equal to f for a = 1. In Fig. [12] 
we use the numerical value of f (|S"5|) and the relation ([54]) 
that connects r' and f to plot a curve in Fig. [T^l which 
agrees with the numerical values of r' for various values 
of a in the same figure. These numerical values of r' 
were obtained by the same numerical methods presented 
in Sec. IIIII and refined by Richardson's extrapolation. 

Because the temperature-dependent parameter a is 
proportional to a density, it should approach its critical 
value a c = 1/3, as |T — T c l^, according to the predic- 
tions of mean-field theory [23, p. 251]. Because surface 
tension vanishes as \T — T r \ 3 ^ 2 in mean- field theory, this 



explains the factor \a— 1 /3 1 3 in ([231) and 



Moreover, 



the results of Varea and Robledo [45j in the mean-field 
approximation show that the ratio of critical exponents 
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FIG. 12. t' as a function of a. The little circles indicate the 
refined r' results calculated directly from various a values. 
The curve is a plot of r' based on a refined calculation of f 
as shown in Table HU 



of line tension and surface tension is 2/3, consistent with 
the ratio of r' in and a' in (|3"5|) . namely 



r 

TP 



oc 



I|2 
3 I 



IIS' 
3 



(55) 



The authors are grateful to one of the reviewers for point- 
ing out this observation. Thus, when the system ap- 
proaches a homogeneous solution, the line tension van- 
ishes more slowly than the interfacial tension. 

We note that a somewhat more general potential, 
namely 



f*(X 1 ,X 2 ,X 3 ) =Y J i^-a i ) 2 {X l -b i ) 2 



(56) 



i=i 



containing the six constants < < 1 and < bi < 1, 
can be mapped onto the potential / in (jig)) . In this 
case, the minima are located at the bulk phases a — 
(ai,b 2 ,b 3 ), (3 = (61,02,63), and 7 = (6i,6 2 ,a 3 ). The 
condition Y^=i Xi = 1 leads to the three constraints 
ai +b 2 + 63 = 1, 61 + a 2 + 63 = 1 and b\ + b 2 + a 3 = 1. 
Regarding the ai to be independent variables, 

bi = (1/2)[1 + a, - a, - a k ] = (1 - Q)/2 + a u (57) 

where i, j and k are all different and Q = Y^e=i a £- Since 

< Q < 3, we have bi > Oj for < Q < 1 and bi < ai 
for 1 < Q < 3. Any choice of the vector (01,02,03) in 
the positive unit cube will lead to bi < 1 but the require- 
ment < bi restricts (01,02,03) to lie within the posi- 
tive unit cube truncated by a pyramid consisting of three 
planes; the apex of the pyramid is located at (1, 1, 1) and 
the other three vertices are located at (1,0,0), (0,1,0) 
and (0, 0, 1). This truncation only restricts (ai, 02, o 3 ) if 

1 < Q < 3. It turns out that the three phases a, /3, 7 
are located at the vertices of equilateral triangles that 
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(a) cii > bi 



(b) a,i < bi 



FIG. 13. Examples of the location of bulk phases (a, /3, 7) 
for the potential /*. (a) For allowed ai and 1 < Q < 3, we 



have ai > bi and the phases are at the vertices of eq uilate ral 
triangles that are magnifications of the Gibbs triangle, (b) For 
< Q < 1, the bulk phases are at the vertices of equilateral 
triangles that are inverted with respect to the Gibbs triangle. 
If Q = 1, we have ai = bi and the triangles degenerate to 
points (bulk criticality) . 



lie within or on the Gibbs triangle and whose sides are 
parallel to the sides of the Gibbs triangle, as depicted in 
Figure fT51 For allowed (01,02,03) and 1 < Q < 3, the 
phases a, /3, 7 are located at the vertices of triangles that 
are magnifications of the Gibbs triangle, as depicted in 
For < Q < 1, the phases lie at the vertices 



Figure 13(a) 



of equilateral triangles that are inverted with respect to 
the Gibbs triangle, as depicted in Figure [13(b) | 

By defining the new variables Zi — 2(Xi — a^)/ (1 — Q), 
which satisfy 5D i=1 Zi = 1, the potential /* becomes 



/* = 



1-Q 



4 3 



(58) 



which has the same form as / in ([48")) . Thus, the potential 
/* is actually a shifted and scaled version of the potential 
/, resulting in r' = t/{B£ 2 ) = [(1 - Q)/2] 2 f . The phases 
merge (bulk criticality) whenever Q = 1. 



V. LINE ADSORPTION 

The Gibbs adsorption equation relates the change of 
interfacial tension to the change of field variables, the 
coefficients being surface adsorptions. In our case, the 
Gibbs adsorption equation for the cv/3-interface is 



da, 



a/3 



TfdM 2 



r^dT 



(59) 



where Tj is the adsorption (surface excess per unit area) 
of the chemical constituent i, and Tt is the adsorption 
related to entropy. Notice that each of the surface ad- 
sorptions, T"^, r 2 1 ' 3 , and T^P , depends on the choice of 
dividing surface, but da a p in the Gibbs adsorption equa- 
tion (|59|) is independent of this choice. Independence of 
location of the dividing surface is one of the properties 
of the Gibbs adsorption equation. A similar relation also 
works for the f3j- and 7cv-interfaces. 
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The Gibbs adsorption equation for diffuse interface 
models have been studied extensively (See Rowlinson and 
Widom [U p 37-38]). As an extension of the Gibbs ad- 
sorption equation, Djikaev and Widom [34j introduce a 
line adsorption equation, which depends on the choice of 
the position of the contact line r. The line adsorption 
equation of Djikaev and Widom is 

c+1 

At = - y^ Aj(f)dfj,i 
i=i 

- (e Q ^dcr a/ 3 -I- epydap*, + e^ a da ia ) ■ (f- r ) , 

(60) 

where A c+ i is the line adsorption corresponding to the 
entropy conjugate to T . This shows that an infinitesimal 
change of the line tension g?t, for a c-component system, 
comes from two parts. The first part is analogous to 
the Gibbs adsorption equation, which is a linear combi- 
nation of the infinitesimal changes of the field variables 
Hi multiplied by the line adsorption A;, that depends on 
the position of the three-phase contact line. The second 
part is the inner product of the difference between r and 
r*o, a specific choice of r, and the summation of the in- 
finitesimal changes of the three interfacial tensions dak 
multiplied by the corresponding unit vector effc along the 
interface k and perpendicular to the contact line. For a 
given value of r*o, the line adsorption equation (|60[) does 
not depend on r. The second term arises because the 
interfaces can change angles as the [ii change. As shown 
by [HI HH, HU , the second term can be eliminated if f~o 
is chosen to lie along a special line. In that case, the line 
adsorption equation (|60[) becomes 

c+1 

dr= -^A,(rOd Ml . (61) 

i=i 

In our case of a symmetric three-phase contact line, 
the second part of the line adsorption equation is 
zero, since the three interfacial tensions remain equal as 
a changes (See the form (f38|) of surface tension) and the 
summation of the three unit vectors is zero. For our case, 
the line adsorption equation (|60|) becomes 

dr = - (A x (r)dMi + A 2 (r)dM 2 + A T (r)dT) , (62) 

where A^ is the line adsorption corresponding to chemical 
constituent i and At is the line adsorption corresponding 
to the entropy. From the form (l62|) of the line adsorption 
equation, it appears that r depends on three variables 
Mi, Mb, and T. However, if we consider the two Clapey- 
ron equations for this three-phase system, 

(p? - p?)dMi + (p% - (4)dM a + (s a - s?)dT = 

{pi - pDdM, + (j4 - pl)dM % + - s^dT = 0, 

(63) 

there is only one independent variable for r, which could 
be Mi, M 2 , or T. For instance, if dr only depends on T, 



we have 

dr = -A c T ff dT, (64) 

where Ajf^ is a linear combination of all Aj and is in- 
variant. By locating our contact line at the center of our 
triangular domain, the symmetry of our potential leads 
to Ai = A2 = 0, so A T ^ = At- A similar simplifica- 
tion also applies for the Gibbs adsorption equation (|5§j) . 
From the Clapeyron equations (f63|) , dMi, dM 2 , and dT 
are linearly related and since (with B£ 2 ^constant) there 
is only one variable a in the problem, we can write 

dr = -A^da, (65) 

where A e Jf = A T dT/da is an effective line adsorption 
corresponding to a. From the relation of r' and f (|54|). 
we calculate 

A *"=£= 2 ^(l) 2 H) f (66) 



VI. SUMMARY AND CONCLUSIONS 

A three-phase contact line in a three-phase fluid sys- 
tem is studied by a mean-field density functional model, 
in which classical sharp fluid-fluid interfaces are replaced 
by diffuse interfaces. The geometry of the system is cho- 
sen to be a prism, where each of its lateral faces is perpen- 
dicular to one of the interfaces and both the cap and bot- 
tom are Neumann triangles. To define a tractable model, 
we assume that the intermolecular forces are short range 
and can be modeled by local densities. The dimension 
of the system is large compared to the interfacial width. 
The excess grand potential of the system is modeled by 
a functional consisting of a highly symmetric three-well 
potential and a gradient energy, which is linear in the 
squared gradients of the three compositions (in terms 
of mole fractions). We assume for simplicity that the 
molar volume is a constant, so there are only two inde- 
pendent densities. We use a variational approach to find 
the governing coupled Euler-Lagrange equations. In the 
far-field limit, where the distance from the contact line 
is large compared to the interfacial width, the transition 
between two bulk phases having different chemical com- 
positions is essentially one-dimensional. Analytically, a 
far-field asymptotic solution is obtained and is used to 
calculate the interfacial tensions. This connects our phe- 
nomenological model to interfacial tensions and to the 
equilibrium angles for classical sharp interfaces. 

Because of the nonlinearity of our Euler-Lagrange 
equations, we cannot find a near-field asymptotic solu- 
tion. Instead, we perform a numerical analysis for a 
symmetric three-phase contact line. By applying a tri- 
angular grid that fills the entire domain, we implement 
a consistent discretization to obtain the discrete Euler- 
Lagrange Equations from the variation of a discretized 
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excess grand potential. To solve the system of these cou- 
pled algebraic equations for the entire domain, we apply 
a successive over-relaxation method and use the asymp- 
totic far-field solutions as the boundary conditions. The 
calculated isoconcentrates (constant mole fractions) bend 
and the effective interfacial width increases slightly near 
the contact line. Close to the outer boundary, the nearly 
parallel isoconcentrates along the diffuse interfaces show 
that our domain size is close to the asymptotic regime, 
so the boundary conditions are sufficient. 

We study the line tension associated with a symmet- 
ric three-phase contact line based on our mean-field den- 
sity functional model, which is the excess grand potential 
over the entire domain diminished by the energies of the 
surfaces extrapolated from the interfacial tensions in the 
far-field. By using the Kerins-Boiteux formula, which 
formulates the expression of the line tension into a sin- 
gle integral, we calculate the line tension and analyze 
the corresponding integrand. Our results show that the 
numerical values of the line tension require a correction 
proportional to the domain size and to the square of the 
grid spacing. To refine our result, we eliminate approxi- 
mately the error associated with the domain size by lin- 
ear extrapolation of the values of line tension from finite 
sizes to zero. Furthermore, we use Richardson's method 
to reduce the error associated with the square of the grid 
spacing to obtain the next level of refinement. The cal- 
culation of line tension based on our mean-field density 
model shows that the value of line tension is negative and 
proportional (a — 1/3) 2 , where a is a parameter in our 
model. We introduce a scaling method to resolve this 
relation for our model. The line tension is proportional 
to (a — 1/3) 2 multiplied by an integral (negative and in- 
dependent of a), in agreement with independent calcu- 
lations for various values of a. In contrast, the far-field 
interfacial tension is proportional to (a — 1/3) 3 . When 
a = 1/3, both the line tension and interfacial tension 
vanish. Physically, this means that the three chemical 
constituents share the same value of mole fraction (1/3) 
and are equally and uniformly distributed over the entire 
system, a single phase. In effect, the interfacial width, 
which is reciprocal to (a — 1/3), is infinite. On the other 
hand, we can either say that the contact line and three 



interfaces vanish or occupy the entire domain. However, 
when a approaches 1/3, the interfacial tensions decay 
faster than the line tension. 

Finally, we relate the change of line tension to the line 
adsorptions by 34]. Thermodynamically, we show that 
there is only one independent field which could be chosen 
as temperature. We are able to link it to the line adsorp- 
tion corresponding to a, since a is the only variable in our 
model (if other coefficients B and t are treated as con- 
stants). Consequently, we find an analytical expression 
of the line adsorption corresponding to a. 

In order to link our model to realistic systems, we 
make the following numerical estimates. For T ~ 300 K, 
typical values of interfacial tension are a few times of 
10 -2 N/m, and interfacial widths are a few A [13. We 
assume a = 1, the interfacial tension a ~ 5 x 10 -2 N/m, 
and the characteristic length I ~ 1 A, corresponding to 
an interfacial width of 2 ~ 3 A. Inserting the form 
for interfacial tension into the form (j54")l for line tension 
yields 



r 
a 



V2 
*"3 



■ft. 



(67) 



We obtain r 0.3 x 10" 11 N. The magnitude of r is at 

the lower end of typical experimental values, which are in 

the range 10" 11 to 10" 9 N k e Jf 0.91 x 10~ n N. 

A crude estimate gives A T r/T ~ 10" 14 N/K. The 

units of At are entropy per unit length. 
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